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ABSTRACT 

A  bubble  transport  equation  has  been  developed  that 
the  bubble  distribution,  as  a  function  of  position, 
velocity,  and  size,  must  satisfy.   An  analytical  model  for 
bubble  transport  in  the  upper  ocean  is  chosen  and  solutions 
are  developed  for  this  model.   For  a  surface  source  only, 
these  solutions  compare  favorably  to  experimental  data  in 
the  near-surface  region.   The  depth  of  this  region  of 
agreement  depends  upon  the  circulation  field  chosen,  but 
is  of  the  order  of  3  meters. 
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z  depth  (z-0) 

v  bubble  velocity 
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<{>  bubble  distribution 
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I.   INTRODUCTION 

Studies  concerning  the  propagation  of  sound  in  the 
ocean  and  under  water  sound  scattering  from  the  sea  surface 
have  shown  a  discrepancy  between  existing  theories  and  ex- 
perimental results.   Because  air  has  a  markedly  different 
density  and  compressibility  than  sea  water  and  because  of 
the  resonant  characteristics  of  bubbles,  the  suspended  air 
bubbles  in  sea  water  have  a  profound  effect  upon  underwater 
sound  and  can  explain  part  of  the  discrepancies  that  have 
been  encountered  (Ref.  1  and  2).   There  are  theories  that 
predict  the  effects  of  air  bubbles  on  absorption,  scatter- 
ing and  attenuation  of  sound  as  well  as  their  effects  on 
sound  velocity  (Ref.  3),  but  there  is  no  theory  at  present 
which  predicts  the  distribution  of  air  bubbles  in  the 
ocean.   Such  a  theory  would  enhance  existing  theories  and 
result  in  better  predictions  on  the  behavior  of  sound  prop- 
agation in  the  upper  ocean. 

A  few  experimental  studies  have  been  conducted  on  the 
distribution  of  bubbles  in  the  ocean.   Blanchard  and  Wood- 
cock were  pioneers  in  this  field,  using  a  bottle  to  scoop 
up  water  from  breaking  waves  and  counting  the  number  of 
bubbles  of  different  sizes  (Ref.  4).   Later  studies  con- 
ducted at  the  U.S.  Naval  Postgraduate  School,  Monterey,  used 
acoustic  techniques  to  determine  the  bubble  population, 
taking  advantage  of  the  fact  resonant  acoustical  scattering 


and  absorption  cross  section  of  bubbles  are  of  the  order 
of  a  thousand  times  the  geometrical  cross  section.   These 
experiments  indicated  that  bubbles  of  radius  60  microns  or 
greater  have  a  distribution  depth  dependence  of  Z  2, 

whereas  smaller  bubbles  follow  closer  the  exponental  law 

-  z/L 
e   '  ,  where  L  is  between  5  and  9  meters  (Ref.  5). 

In  this  paper,  development  of  a  bubble  transport 

equation  has  been  outlined  that  the  bubble  distribution,  as 

a  function  of  position,  velocity  and  size,  must  satisfy. 

An  analytical  model  for  bubble  transport  in  the  upper  ocean 

is  chosen  and  solutions  are  developed  for  this  model. 


II.   TRANSPORT  EQUATION 

A  heuristic  approach  is  used  here  to  develop  the 
bubble  transport  equation,  since  the  results  are  the  same 
as  with  a  more  rigorous  approach  and  the  mathematics  is 
less  involved  (Ref.  12).   Consider  the  distribution 

y  (r ,v ,1 , t)d°rd  vdl  which  defines  the  number  of  bubbles 

•z 

with  radii  in  dl  about  1  and  with  velocity  in  d  v  about  v 

which  are  in  a  volume  d  r  about  r  at  time  t.   If  collisions 
among  the  bubbles  are  neglected,  then  the  following  equa- 
tions can  be  written  for  the  seven-dimensional  space: 

|X(r,v,l,t)d3rd3vdl  =  [net  flux  in  +  S(f,v,l,t)J  d3rd3vdl  (1) 
where  S  (r ,v , 1 , t)d  rd  vdl  is  a  distributed  source  of  bubbles 
introduced  in  d^r  about  f  at  time  t  with  radii  in  dl  about 
1  and  velocity  in  d  v  about  v. 

Referring  to  figure  1,  which  represents  a  volume  ele- 
ment in  real  space,  the  net  flux  in,  per  unit  volume,  for 
the  x-direction,  is: 

fluxin-fluxout  _vx  (x)  V  (x)  dydz - vx  (x+dx)  V  (x+dx)  dydz   3  (vyy) 
volume  dxdydz  ="8x 

Similiarly,  the  net  flux  in  per  unit  volume  for  the  y  and 

z-directions  is  -  —7 LJ-  and  -  — ^ }-t    respectively. 

3y  dz 

Combining  the  above  results,  the  total  net  flux  in  per  unit 
volume  for  real  space  is  -  Vr,(vvF),  where  the  subscript  r 
on  the  del-operator  means  the  operation  is  taking  place  in 
the  real  space.   A  similiar  approach  is  used  to  obtain  the 
total  net  flux  in  for  velocity  and  size  space:  -  V  • (a¥) 


vx(x)¥(x)dxdz 


(Flux  in) 


vx(x+dx)y(x+dx) 
dydz 

>(Flux  out) 


Figure  1 

Elemental  Volume  Showing  "Bubble" 
Streaming  in  Real  Space  (x,y,z). 


for  velocity  space  and  -  8 (R4Q  for  size  space,  where  R=dl/dt 

The  total  net  flux  in  for  the  seven-dimensional  space  under 
consideration  is  the  sum  of  the  net  flux  in  for  real, 
velocity,  and  size  spaces.   Substituting  this  into  equation 
(1)  gives  the  general  bubble  transport  equation: 

11   =  _  V  M)    -   V  •  (a¥)  _  3(RY) 

at     r       v      91 —  +  b     (2) 

The  remainder  of  this  paper  will  consider  the  special 

case  of  steady-state  with  no  distributed  sources,  for  which 

equation  (2)  can  be  written 

v.VH'  +  a.V  <F  +  Ml  =  -  fv  .a  +  Ml  y      C3") 
r       v     81     L  v     9i  J        K    J 

since  f  and  v  are  independent  variables.   Although 
distributed  sources  will  be  neglected  in  this  thesis, 
surface  sources  can  be  introduced  as  boundary  conditions. 
It  should  be  pointed  out  that  the  non-zero  right  hand 
side  of  equation  (3)  is  peculiar  to  the  bubble  transport 
problem  and  arises  becouse  a  is  a  function  of  v  and  R 
is  a  function  of  1. 


III.   ANALYTICAL  MODEL  FOR  TRANSPORT  EQUATION, 

Before  solutions  to  equation  (3)  can  be  obtained,  a 
realistic  model  for  the  bubble  transport  equation  must  be 
developed,  with  particular  attention  to  formulating  ex- 
pressions for  the  bubble  acceleration  (a),  the  time  rate  of 
change  of  the  bubble  radius  (R)  and  the  circulation  field 
(V). 

A.   MODEL  ASSUMPTIONS 

The  model  used  in  this  development  assumes  Reynold's 
numbers  less  than  one,  which  is  the  Stokes'  Law  Regime. 
With  regard  to  this  assumption  two  comments  should  be  made. 
First,  for  small  Reynold's  numbers,  bubbles  will  behave 
like  solid  spheres  with  little  or  no  distortion  and  secondly, 
this  places  an  upper  limit  of  about  150  microns  on  the  size 
of  bubbles  that  can  be  considered.   For  Reynold's  numbers 
much  greater  than  one,  (larger  bubbles)  the  resistance  of 
the  turbulent  wake  will  have  to  be  included  in  the  equation 
of  motion,  along  with  Stokes'  drag  force  (Ref.  6). 

Surface  tension  has  been  neglected  in  the  model,  so  the 
pressure  inside  the  bubble  is  equal  to  the  hydrostatic 
fluid  pressure.   Since  the  pressure  from  surface  tension 
2y/l|  where  y    is  approximately  74x10"^  NT/m  for  sea  water, 
this  assumption  will  introduce  an  error  of  less  than  5 
percent  for  radii  greater  than  30  microns.   Thus,  where  the 
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Reynold's  number  places  an  upper  bound  on  the  radius,  neg- 
lecting surface  tension  places  a  lower  bound  on  the  bubble 
radius  under  consideration. 

The  volume  of  a  gas  bubble  in  a  liquid  varies  primarily 
through  the  effect  of  two  influences.   First,  gas  contin- 
uously leaks  out  of  the  bubble  to  dissolve  in  the  surround- 
ing liquid.   If  the  bubble  was  stationary,  the  ensuing  loss 
of  mass  would  necessarily  lead  to  its  ultimate  disappearance 
On  the  other  hand,  as  a  bubble  assends  or  desends  in  the 
liquid,  the  ambient  hydrostatic  pressure  decreases  or  in- 
creases and  so  does  the  pressure  within  the  gas  bubble,  thus 
effecting  the  size. 

The  effects  of  the  interaction  between  gas  diffusion 
and  change  in  hydrostatic  pressure  can  best  be  seen  by 
looking  at  the  ideal  gas  law: 

^■n-1-5  =  rTN/P 

where  N  is  the  total  number  of  gas  molecules,  r  is  the  gas 
constant,  T  is  the  temperature,  and  P  is  the  pressure.   Now, 
taking  the  time  derivative  of  the  ideal  gas  law  equation 
and  assuming  isothermal  changes, 


4tt12r  =  IT   pdN/dt  -  Ndp/dt 
P 


and  rearranging 


1    dN/dt         dP/dt 

i      N  p 
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The  above  development  assumes  a  single  gas  inside  the 
bubble  vice  some  mixture  of  gases.   In  equation  (4),  the 
first  term  is  a  result  of  gas  diffusion  whereas  the  second 
term  brings  in  the  effects  of  compression. 

Due  to  the  relative  motion  between  the  water  and  bubble, 
assume  that  the  composition  of  the  water  in  respect  to  dis- 
solved gas  is  everywhere  uniform  and  saturated  at  atmospheric 
pressure,  P  ,  except  in  a  thin  shell  surrounding  the  bubble. 
Assume  that  the  gas  in  the  bubble  is  uniform  right  up  to 
the  gas-water  interface,  and  that  the  liquid  in  contact 
with  the  bubble  is  saturated  at  the  pressure,  P,  inside  the 
bubble.   Further  assume  that  diffusion  gradient  is  uniform 
throughout  the  thin  shell.   With  these  assumptions,  Fick's 
First  Law  can  be  used,  which  relates  the  amount  of  material 
which  diffuses  through  an  unit  area  to  the  product  of  the 
concentration  gradient  and  diffusion  constant: 


dWdt    =  -  64irl2(P-Po) 


In  this  expression  6  is  a  proportionality  constant  which 
depends  on  d,  the  thickness  of  the  shell;  on  A,  the  dif- 
fusion constant  of  the  gas  in  the  water;  and  on  a,  the  solu- 
bility of  the  gas  in  the  water,,  in  accordance  with  the 
relation : 


6  =  Aa/d 
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Using  Fick's  First  Law  with  the  hydrostatic  pressure  rela- 
tionship, P=P   -  Pwgz,  equation  (4)  becomes: 

v 

r  =  6rT  — - +  -  — - (4) 

(cq-z)     3  (a1-z) 

where  a^=  P0/pwg  and  vz  is  the  bubble  velocity  in  the 
z-direction.   In  reference  7,  the  expression  6rT  was  deter- 
mined experimentally  to  be  approximately  a  constant, 
10'6  m/sec. 

B.   EQUATION  OF  MOTION 

When  air  bubbles  are  entrained  in  a  circulation  field, 
there  are  three  forces  acting  on  them  along  with  the  bouyant 
force.   Figure  2  illustrates  an  air  bubble  in  a  circulation 
field  and  shows  graphicly  the  various  forces  which  need 
to  be  considered  in  order  to  determine  the  bubble 
acceleration . 

If  the  bubbles  were  completely  entrained  by  the  fluid, 
i.e.,  bubble  velocity  =  circulation  field  velocity,  they 
would  experience  the  same  force  as  would  the  fluid  enclosed 
in  the  same  volume,  pwVdV/dt.   However,  complete  entrain- 
ment  can  only  occur  if  p  =  pw ,  which  is  far  from  being  true 
for  air  bubbles.   Since  the  bubbles  are  only  partially  en- 
trained, the  fluid  flows  past  the  bubbles  creating  a  rela- 
tive velocity  (u)  and  thus  introduces  a  drag  force.   As 
mentioned  earlier,  small  Reynold's  numbers  have  been  assumed 
and  thus  the  drag  acting  on  the  bubble  can  be  described  by 
Stokes'  Law, 
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*x 


Bouyant  Force 

A 


Entrainment  Force 


Drag  Force 


^Adhering  Mass 
Force 


Figure  2 

Free  Body  Diagram  of  a  Bubble, 

Showing  the  Forces  Acting  on 

the  Bubble  in  a  Circulation  Field 
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FDrag  =  "  6^lu- 

The  third  force  acting  on  the  bubble  is  the  innertial 
adhering  mass  force.   When  a  spherical  shape  moves  through 
the  water,  it  drags  along  an  amount  of  water  equal  to  a 
fractional  amount  of  the  volume  of  the  sphere.   In  the 
equations  that  follows,  sigma  (a)  represents  this  fractional 
part,  and  in  reference  (8),  sigma  has  been  derived  as  0.5 
for  a  sphere. 

Referring  to  figure  2  and  using  Newton's  Second  Law  of 
Motion,  the  equation  of  motion  for  an  air  bubble  in  a  cir- 
culation field  can  be  written  as: 

•  ■  pVg  -  PWV||  -  crVpw||  -  6,mlu  ♦  (pw-P)ffgk. 

Using    P<<PW   and    dV/dt    =    (v-Vr)V   yields 

a   =    -    °L_  a   +   £k   +   ilPlCv-VrlV  (5) 

2  o  o 

where  a  =  9/  2  n/  crp  and  u=v-V. 

W 

C.   CIRCULATION  FIELD 

Little  is  known  about  the  circulation  field  in  the 
upper  ocean.   Although  some  work  is  being  conducted  at  pre- 
sent at  the  U.S.  Naval  Postgraduate  School,  Monterey.   In 
reference  (6) ,  Levich  develops  a  two-dimensional  equation 
which  describes  the  particle  motion  by  wave  action  in  an 
ideal  liquid,  and  is  the  basic  form  used  here  to  describe 
circulation  field: 
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V  =  VQe+  zsin  kxi  -  V  e+kzcos  kxk.  (6) 


This  equation  gives  a  periodic  form  for  the  horizontal 
direction  as  well  as  an  exponental  decay  in  depth,  which 
might  be  expected  for  the  physical  situation  in  the  upper 
ocean.   It  is  also  well  to  note  that  the  condition  for  an 
imcompressible  fluid  has  been  satisfied  with  this  equation, 
i.e. ,  V«V=0. 

Figure  3  shows  the  cell  structure  implied  by  this  cir- 
culation field.   For  this  model,  any  two  adjacent  cells  are 
identical  to  every  other  pair  of  adjacient  cells  that  are 
separated  from  the  first  pair  by  an  integer  number  of 
wavelengths.   This  implies  that  the  boundary  conditions  must 
be  periodic  unless  unphysical  surface  sources  are  to  be 
assumed . 
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kx-> 


z=o 


,  2ir 

K  = 

A 


Figure  3 

Cell  Structure  of  the  Circulation 
Field  Model 
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IV.  SOLUTION  TO  TRANSPORT  EQUATION 

In  this  section,  a  look  is  taken  at  the  bubble  transport 
equation  with  the  model  that  was  presented  in  the  previous 
section.   The  main  interest  will  focus  at  determining  the 
bubble  distribution  as  a  function  of  depth  and  bubble  size. 

Substituting  into  the  bubble  transport  equation  the 
expressions  for  a  and  R  from  the  preceeding  section  results 
in  the  following  equation: 


v-VtY  +  a-VvY 


v 


+  R^/31  -  [3°7l2  "  ^/SCc^-z)] 


* 


By  defining 


<J>(z,x,vz,vx,l)  = 


00 

t    f 


00 


y(r,v,l)dydv 


y 


and  using  the  two-dimensional  circulation  field,  the 
above  equation  can  be  integrated  over  the  y  and  v   directions 
and  the  problem  reduces  to  solving  the  partial  differential 
equation : 


vxli  +  vz!4  +  axil  +  a  dj     Rdl  m 
^x      9z      3VX     3vz+  TT 


o      V 

2a    z 
2  "  3Xa1-z) 


with  the  boundary  condition  <j)(0,x,v  ,v   1)  =  <f>   (0>xo>vzoj 
vXQ , 1  )  corresponding  to  a  surface  source. 

Equation  (7)  is  a  first  order,  linear  partial  differen- 
tial equation  in  5  variables  and  can  be  solved  by  the  method 
of  characteristics  (Ref.  9).   Using  this  method,  the  partial 
differential  equation  is  equivalent  to  the  following  5  char- 
acteristic ordinary  differential  equations  in  terms  of  the 
parameter  z,  which  was  chosen  for  mathematical  convience: 
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Characteristic  Diff.  Eqns 

Z<0 

dx   *x 

3Y  '  vz 


dv. 


37*= 


Characteristics 
z  =  0 
x=x„ 


vx=vx0 


dv7=  a^ 

dz    v„ 


vz=vzo 


dl 

dz 


R 


1  =  1. 


d£ 
dz 


2a 


1  vz     3(ai-z)J 


<J>     A  =  <j)n(x  ,v    ,v    ,1  ) 
y  to  ^  0  »  xo  '  zo  '  oJ 


A  close  look  at  the  characteristic  ordinary  differential 
equations  reveals  that  they  are  simply  the  bubble  dynamic 
equations . 

Solving  these  equations  for  <$>   would  be  straight  forward, 
except  the  equation  for  v   is  non-linear  and  analytic  solu- 
tions have  been  found  only  for  very  special  cases  (Ref.  10). 
Thus,  approximate  methods  must  be  considered  to  obtain  v  . 

The  method  used  for  this  paper  to  solve  for  v   was  to 

assume  v   =  V   +  u  .   Appendix  A  outlines  this  method  with 
z     z     z     rr 

the  solution  for  uz  as: 


uz  =  e 


"Jo  I7^ 


dz 


fz1 

a 

0   2 
1  v, 


dz 


av_    alv7  8x  dz/ 


dz  '  +u 


zo 


Since  the  drag  acceleration  on  the  bubble  is  large,  the 
bubbles  will  be  entrained  by  the  fluid  within  the  first  few 
centimeters  of  travel  and  u   thus  approaches  a  terminal 
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velocity.   Referring  to  the  computer  output  for  u  ,  shows 
that  the  assumption  V   >>  u   can  be  made  and  also  u   is 

Li  £->  L, 

independent  of  V   ,  thus  implying  that  v   is  independent 

zo 

By  a  similiar  approach,  it  is  shown  in  Appendix  A  that 

v  z   Vx .   Now,  the  characteristic  equation  involving  x  can 

be  written  as: 

A  V 

dx  ~    vx 

a?     vT  ' 

and  using  the  chosen  circulation  field, 

x  *  1/k  sin   (e~*zsinkx0) . 
This  equation  implies  that  the  bubbles  follow  trajectories 
in  the  fluid  and  form  streamlines  as  was  indicated  on  figure 
3.   Thus  choosing  a  x   one  can  describe  the  path  which  the 
bubbles  follow. 

The  solution  of  the  characteristic  equation  for  1  is 
straight  forward  with 


1/3 


fz 


2/3 


l--U-z/a1)-A'J|i0  +  iH  s'Cl-z'/cn)    dz 

x  L      a  J  o     ^z 

The  second  term  in  this  equation  is  a  result  of  gas  diffu 
sion  and  limits  the  bubble  lifetime  and  hence  the  depth 
that  surface  generated  bubbles  can  penetrate. 

From  the  last  characteristic  equation,  <J)  takes  the 
general  form: 

*(z»xo»vxo»vzo»1o)  =  (^0^X0>VXO'VZO'10^ 


(l-z/a1)1/3  EXP 


fZ2a   dz' 
°l2vz   J 
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where  <J>  as  a  function  of  (x,z,v  v7,l)  is  obtained  from  the 
relationships  between  the  characteristics  and  their  variables 
Assume  a  surface  source  cj>   that  is  separable, 

4>  (x  ,v   ,v   ,1  )  =  X(x  )Y  (vYn)Y(v   )L(1  ), 
o   o   xo   zo   o       o   x   xo'  v  zo;  l  oJ  > 

and  for  convience,  let 

Yx(vxo)  ■  «CvX0-vx)   !   Vv20)  "  «fvzo^z'- 
Since  v   and  vz  are  nearly  independent  of  vxo  and  v_0,  this 
choice  will  not  strongly  affect  the  solution.   Therefore, 


<j>(x,z,v  v  ,1)  =  X  (x  )L(1  )<S(v   -v  )6(v   -v  ) 

X   Z         0   0      O     XO   X     zo  1 


1/3 


2a   dzt 
0l2v 


Now,  define  the  "bubble  density"  as  a  function  of  position 
and  size, 


t°P 


to 


<\> (x ,z ,v  ,v  ,l)dv  dv 
x   z      X   z 


$(X,Z,1)   E 

and  integrate  the  expression  for  <$>  over  vx  and  v  : 


*(x,z,l)  = 


X(x  )L(1  )6(v   -v  )6(v   -v  ) 


XO   X' 


ZO    Z    L_ 


(1-z/ 


aj) 


/3 


EXP 


fz2a 
ol^v, 


dz« 


dvv  dv 
x    z  > 


To  take  advantage  of  the  delta-functions,  the  variables  of 
integration  are  changed,  dv  dv   =  Jdvxo  dv   ,  where  J  is 


the  Jocobian: 


J  E  3CVX»Vz) 

^vxo'vzo^ 
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The  equation  of  interest  now  has  the  form: 

$(x,z,l)  =  X(x  )L(1  )(l-z/a1)1/3  EXP 
Because  of  the  number  of  mathematical  steps  involved  in 


2a   dz'J   (8) 
oi2Vz 


deriving  the  Jocobian,  the  development  is  found  in  Appendix 
B.   The  final  approximate  result  is: 


J 


v 


zo 


rZ 


V. 


2a 


<>i2v. 


dz1  EXP 


J«-ea«- 


where , 


=  +12_  (v  -V.)|ft 

l^vz    z   z  ai 


-Z'  (l-z'/aj 


dz 


O   V- 


Thus,  the  equation  for  the  distribution  of  bubbles  with  one 
initial  surface  velocity  can  be  written  as  (Ref.  12): 

Hx,z,l)  *  X(x0)L(l0)(l-z/al)1/3^  ej^dz'     (9) 
Looking  at  equation  (9),  the  term  eQdz'  will  be  the  pre- 
dominate feature  in  solving  for  <S> .   Two  things  should  be 
noted  about  £,   First,  £  is  seen  to  be  proportional  to  1/V  3 

when  the  approximation  v  ~V   is  used.   This  result  illustrates 
rr  z   z 

the  importance  of  selecting  a  good  circulation  field  model. 
Secondly,  £a6RT  and  shows  that  gas  diffusions  is  very 
important  in  the  solution  for  the  bubble  distribution. 
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V.   RESULTS 

A  computer  program  was  written  for  calculating  equation 
(9),  to  determine  if  the  results  were  physically  realistic. 
In  this  program,  the  assumptions  were  made  that  X(x  )=1  and 
L(1Q)=K1  "7/2  (Ref.  11).   For  the  latter  assumption,  K=l 
was  used  in  the  computions  since  the  interest  was  in  the 
shape  of  the  distribution  curve  rather  than  the  actual  den- 
sities.  Calculations  for  three  circulation  field  veloci- 
ties were  made  and  the  following  graphs  have  been  plotted: 

Zmax  vs-  1o 

$  (0,z,l)   vs.   Z   for  various  radii 

$  (0,z,l)   vs.   1   for  various  depths. 

Z     is  the  maximum  penetration  depth  for  a  given  size 

Jit  ci  A. 

surface  generated  bubble.   All  calculations  were  made  for 

x  =0,  which  implies  a  straight  line  trajectory  downward. 

If  some  other  x   had  been  chosen,  line  integration  would 

o  '  & 

had  to  be  performed  on  the  characteristic  curves  in  5-space 

(x,z>vx,v3,l) . 

The  numerical  technique  used  in  computing  the  integrals 
was  10-point  Gaussian  Quadrature,  which  approximates  the 
integrand  by  a  Legendre  polynomial  series.   Although  this 
technique  was  fairly  accurate  for  most  of  the  integrals,  it 
could  not  be  used  in  computing  the  relative  velocity  (uz) . 
For  this  calculation,  it  is  noted  that  the  integrand  is 
essentially  zero  except  for  the  last  few  percent  of  the 
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interval  (0,z)  where  it  increases  dramatically.   U  was 
therefore  evaluated  by  using  double  precision  32-point 
Gaussian  Quadrature  on  the  latter  portion  of  (0,z). 

In  computing  the  relative  velocity,  the  approximation 


rz 


a 


1  vz 


dz'  *  ^ 


12V0 


-kz     n 
ze    ,  x=0, 


was  used  which  is  accurate  provided  |Z|<  10  meters.   Also, 
since  the  calculations  involved  taking  the  difference  be- 
tween two  such  integrals,  very  little  error  was  introduced. 

It  is  appropriate  now  to  make  a  few  comments  about  the 
results  that  have  been  plotted.   Figure  4  is  a  graph  of 
maximun  depth  versus  inital  bubble  size.   This  graph  shows 
that  bubbles  have  a  finite  life  time  as  a  result  of  gas  dif- 
fusion and  that  smaller  bubbles  introduced  at  the  surface 
will  not  penetrate  as  deep  as  larger  bubbles. 

Figures  5  through  7  are  plots  of  $  (0,z,l)  versus  depth 
for  various  bubble  radii.   These  graphs  show  the  influence 
that  the  circulation  field  velocity  has  on  how  deep  the 
bubbles  will  penetrate.   It  should  also  be  noted  that  even 
though  the  surface  density  falls  off  as  1  "3.5^  there  is  a 
certain  depth  where  the  larger  bubbles  become  more  populous 
than  the  smaller  bubbles  because  of  gas  diffusion. 

Figures  8  through  11  are  plots  of  $  (0,z,l)  versus 
radius  for  various  depths.   The  main  feature  about  these 
graphs  is  that  the  initial  distribution  appears  to  persist 
to  some  depth  which  is  a  function  of  the  circulation  field 
velocity. 
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Figure  4 

Curves  for  Maximun  Penetration  vs.. 

Bubble  Radius  for  Various  Circulation 
Field  Velocities 
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Figure  5 

Bubble  Density  vs  Depth  for  Various 

Size  Bubbles,  v  =-0.5  m/sec. 
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Figure  6 

Bubble  Density  vs  Depth  for  Various 
Size  Bubbles,  v  =-1.0  m/sec. 
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Figure    7 

Bubble  Density  vs.  Depth  for  Various 
Size  Bubbles,  v0=-1.5  m/sec. 
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Figure  8 


Bubble  Density  vs  Bubble  Radius 
for  Various  Depths,  v  =-0.5  m/sec 
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Figure  9 

Bubble  Density  vs  Bubble  Radius 

for  Various  Depths,  v  =-1.0  m/sec 
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Figure  10 

Bubble  Density  vs  Bubble  Radius 
for  Various  Depths,  v0=-1.5  m/sec 
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VI.   CONCLUSIONS 

There  is  little  experimental  data  available  to  compare 
the  results  to,  however,  in  reference  (5)  where  measurements 
were  taken  acoustically,  general  shapes  for  bubble  distri- 
butions are  given.   One  should  note  that  the  experimental 
data  in  reference  (5)  is  not  for  isoradius  bubbles  since 
the  resonant  frequency  is  a  function  of  depth,  however,  a 
good  comparision  can  be  made  within  a  few  meters  of  the 
surface . 

Although  the  bubble  distribution  as  computed  does  not 
follow  a  simple  mathematical  function,  the  results  for  the 
first  few  meters  can  be  approximated  by  exponental  laws, 
depending  on  the  circulation  velocity  and  bubble  size. 

Looking  first  at  figure  5,  the  distribution  drops  off 
rapidly  with  depth  pointing  to  the  fact  that  for  small  cir- 
culation velocities,  a  distributive  source  is  needed  since 
the  depth  penetration  is  small.   For  V  =-1.0  m/sec,  the 
curves  can  be  approximated  for  the  first  3  meters  by  the  ex- 
ponental laws,  e~z'   for  small  bubbles  and  e~z'   for  large 
bubbles.   The  curves  for  V  =-1.5  m/sec  follow  fairly  close 
e"z/9  for  the  smaller  bubbles  when  z^-3.5m  and  e"z/15  for 

larger  bubbles  when  z>-4.5m.   The  curves  for  V  =-1.0  and 
to  o 

V=-1.5m  /sec  compare  favorably  with  the  experimental  data 
mentioned,  especially  for  the  smaller  bubbles. 

Even  though  there  is  fairly  good  agreement  in  the 
distribution  near  the  surface,  a  definite  discrepancy  does 
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occur  at  the  deeper  depths.   Probably,  the  major  reason  for 
this  discrepancy  is  the  fact  that  no  distributive  or  bottom 
sources  were  considered  in  the  calculations.   The  data  in 
reference  (5)  was  taken  in  shallow  water  and  calm  seas,  so 
inclusion  of  such  sources  into  the  model  should  improve  the 
agreement  between  the  solution  and  the  experimental  data. 
Another  possible  reason  for  the  discrepancy  is  that  a  large 
portion  of  the  Jacobian  was  assumed  small  and  neglected. 
This  term  was  not  looked  at  in  detail  and  could  be  of  some 
importance  at  the  greater  depths. 

Since  the  results  do  agree  in  part  with  experimental 
data,  the  approach  taken  in  this  paper  does  appear  feasible 
for  determining  the  bubble  distribution  in  the  upper  ocean. 
For  further  studies  in  this  area,  it  is  suggested  that  solu- 
tions to  the  bubble  transport  equation  be  extended  to  include 
calculations  for  x^O,  models  with  distributed  sources  and 
surface  tension  and  better  models  for  the  circulation  field 
and  gas  diffusion. 
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APPENDIX  A 

Solutions  for  the  Relative  Velocities 

In  this  appendix,  solutions  are  obtained  for  the  rela- 
tive velocities  from  the  characteristic  equations.   Consider 
first  the  characteristic  equation  for  v  ,  which  is  non-linear 
Write  v  =V  +u  ,  differentiate: 


dv 


-,  =  dVz  +  duz 


dz     dz     dz 
Substituting  this  result  into  the  characteristic  equation 
for  vz  gives 


du. 


Hz-  +  i2 


a  uz 


g 


v- 


av- 


1 

\      8Vz_    9Vz 

rz   3Y   +  8Z 


A  formal  solution  to  this  equation  can  be  written: 


u. 


dz 


u 


n—    dz' 

012VZ 


zo 


p 

'z  * 

e-1 

LJ 

0 

fz  »  —7-  0Qz  r 
l2v  2 

o      L 


l(vx    8vz 


3v- 


aVz  +  a\vz  9x   +  9z 


J 


or,  rewritting, 


rz  rz_a 

u7  =  J  e  z.,2 
z    J  0    1  v. 


dz 


11  _. 


av7   alv, 


Sv 


8x 


z  + 


9v 


r  z 


5^  dz'  +  uzoe  J 


a 


o  T  dz 
1  v7 


For  computational  purposes,  the'  approximation  v   =  V  will 

Lt  l-t 

be  made  where  appropriate,  since  the  above  expression  is  an 
implicit  equation.   This  approximation  turns  out  to  be  valid 
as  indicated  in  the  Results  Section. 
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In  order  to  solve  the  characteristic  equation  for  v  , 
the  same  approach  is  used  as  described  above,  with  the 
result : 

r~ 
^—  "^TT  —    dz 


U 


"'oJl2^ 


1  dvx 


1   8V, 


rz  (Z1  r 

— 
0  L.2   "av^8x 
II  v„    <- 


dz' 


vx  av 


v7  8x 


dz 


1 

'  +  u 

xol, 


[_a^z 


A  close  look  at  this  equation  reveals  that  ux  is  approximately 
zero  since  the  exponental  term  shows  a  sharp  decaying 
feature  and  the  circulation  field  chosen  does  not  c1  ige 
rapidly.   Thus  the  approximation  vx=Vx  can  be  used  i   !:he 
calculations . 
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APPENDIX  B 


Development  of  the  Jacobian 


This  section  is  devoted  to  the  mathematical  manipula- 
tions involved  in  developing  an  expression  for  the  Jacobian 
required  for  equation  (8) .   The  following  are  the  functional 
relationships  among  the  variables  and  characteristics  that 
need  to  be  considered: 

^o^zo'V 

V  ,v   ,1  ) 
xo   zo   o 

V  ,v   ,1  ) 
xo'  zo '  0J 

vxo'vzo'V 

vxo  *  vzo  »  o  '  xo-' 

v   ,v   ,1  ,V   ). 
xo   zo   o   zo 


vx    = 

vx 

:z,x0, 

V 

z 

V 

z 

[z,xo, 

1    = 

1 

[z,xQ, 

X 

X 

[z,xe, 

Vx   ■ 

Vx 

Cz,xQ, 

V      = 

V 

'z  ,x    , 

For  mathematical  convience ,  define: 

Jx  E  8vx/8vxo   ;   J2  =    8vx/9vzo   ;   J3  =  3vz/8vxo 

J4  =  dvz/dvzo       ;   J5  =  3x/3vxo   ;   J6  =  3x/8vzo. 

Thus  , 

9(vX>vz) 
=  9(vX0,vZ0)  =  J1J4  "  J2J3- 

Now,  using  the  characteristic  equations  for  v  ,v   and  x 

x   z 

and  noting  that,  for  example: 

3     (&)- 


dv 


xo 


dJ-i 
dz   ' 
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the  following  equations  can  be  written 

r 


dJi 


l2v2   °   Vz 


Ji  + 


a  (vy-Vv)    1+a  k  V. 


.  JC IX- 

vz2    + 


V 


a   v. 


x   J. 


_a_  3VX     Ha  Xx  k  3V£_   1+a  k  dV^   +  2a   (vY-VY)31 

8vxo  "   a   vz    3vxo     a     3vxo   ^v z       3vx0 


l2v. 


dJ 


3z 


2  = 


a 

'7vz 


1+a  kVz 


a  v. 


J. 


_£  (vx~Vx)    1  +  °  v 
I2   V 


:  Vz 


a   v  2 
z 


J4 


+   _a dV^      _    1+a   v^  k   SiV^     +liak8Vx+2a      (vx-Vx)31 


!2vz  ^zo     a   vz    8vzo 


0     9vzoVvz 


3v 


zo 


dJ? 
3z 


a   V 


2   ? 
1  vz2 


"£  2_  i±£Xx  kVx 

avz    a  vz2   X_j 


J   + 
3 


1+a  k  V 


a     v 


*  Ui 


a    3  y7    1+a  Vv    3VX    1+a    3V7    2a  (vz-V z)31 

5   k  »vX0+13Vz 


l2vz  8vxo 


a  vz  k  3vxc 


3v 


xo 


dJ4  = 
dz 


_a  Vz_ 
12  vz2 


g 

avz2 


1+a  vx 

a  vz2 


kV, 


J 


J4  + 


1±£  k-Yxlj. 

a   vz    J  ' 


a    3V7     1+a  vx    3VY     1+a  v  3VZ     2a   (v  -V7)31 
l\z  ^zo    a   vz  k  9vzo     a     3vZo    l3vz       8vzo 


For  ease  o£  manipulation,  define 


A  =  .  _a    _  1+a  kV^  .   B  E  _a_  (vx-Vx)  +  i±£  kVz  v^ 

2      a   v,   '      J   v72      a      v72 
I  vz        vz         1     ^  ^ 


r  =  -  a       S.    1+a  kVxvx     n 
c  -  7T   v  ™_? 

1  vz2  vz 


av„2    a   vz2 


1+a  kV, 


a   v- 


E  = 


a  kVz  _  1+a   vx  k2yx  _  l+o   2 


l^v. 


a   v. 


k"V, 


a  kVY  _  1  +  a   vv   kV,    1  +  a  k2V, 


2,r        a 

1  vz 


v* 

z 
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qi  e  ^a   (vx-Vx)   ;       2a    (v.-V^ 
l3vz  =  l3Vz 

The  six  equations  can  now  be  written  as: 


dJn 

3^  =  AJt  +  BJ3  +  E  J   +  q  £ 

D     i  dv 


xo 


dJ2  m   AJ2  +  BJ4  +  E  J6  +  q   91 

dz  -1  3v 


zo 


dJ- 


81 


g^  =  CJ3  +  DJX  +  F  J5  +  q2  ^ 


xo 


£Ii.  =  CJ4  +  DJ2  +  F  J6  +  <\2   li 

oz       4      ^  av 


zo 


^5.  =  L_  j   -  Vx  j 

dz     vz   1    ^— j   3 


oz    vz   2    vz2  J4 


Taking  the  moments  of  the  first  four  equations 


dJ- 
dT 


J4  «±L 


=  AJlJ4  +  BJ3J4  +  EJ4J5  +  qi3j_  J4  =  ^(j1j4)-jidj4 

3v, 


xo 


J3  dJ^    AJ2J   +  BJ„J4  +  EJ^J6'+  qn81   J3  =  d  ( J2 J3) " J2dJ3 
dz   ~    •        *  ^vzo      dz*         aT 


j   dJ^   CJ2J3  +  DJnJ2  +  FJ2J5  +  q?^l   J   =  d  (J  JJ-J.dJ? 


dz 


9vxo 


o^ 
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Ji  d^4  =  CJlJ4  +  DJltJ2  +  FJ1J6  +  qi|l   J1=  d^J^KJ^ 


3vzo  x      aT 


Now,  combining  the  above  equations  give  the  result: 


£   -  (A+C)J  +  \ 


where 


11 

3v 


J 


xo 


91 


J 

3 
zo   J 


q0  I  31   J     81    Jn 
zo       xo 


Q  =  eCJ4J5  "  J3J6}  +  F(J1J6  "  J2J5} ' 

Using  the  equation  for  1  that  was  developed  earlier  and  taking 
the  partial  derivatives  gives: 


ii 

3v 


6rt 


tz 


xo 


and 
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3v 


a- 


6rt 


'o 


-  2/3 

•'  (l-z'/aQ     J3  dz 

vz2 


(l-z/ai)"1/3 


zo 


a- 


Zz'  (1-z'/oj)"2/3Ja  dz'   (l-z/a-,)"1/3 
o      v  2  J 


If  gas  diffusion  had  been  neglected  in  the  development,  the 
above  two  equations  would  have  been  zero  and  the  solution  to 
the  transport  equation  would  become  unphysical. 

Substituting  the  above  results  into  the  differential 
equation  for  J  results  in: 


dJ   _(A+C-£)J   +    q. 


5rt 


rz 


a 


IJ 


z'  (l-z,/ai)"2/3(J   J    '     "    J    'J    )dz' 


v    2 
z 


3    4 


3      4 


_q. 


rz. 


6rtpz'  (l-z'/aT) 
alJo  v?2 


2/3 


[(JiV    -    J2J3,)-(Jl,J4  "  J2'J3) 


dz 


r 


Q 
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where  now, 


5  =  q 


6rt 
2   on 


r  z 


z«(l-z'/a1)-2/3d; 


v  2 
z 


For  simplification  and  as  a  first  approximation  assume 


HI  =  (A+C-»J 

and  the  solution  to  this  differential  equation  is 


J  =  EXP 


fz 


(A+C-Odz' 


since  J  =  1  at  z=o. 
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